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2-D Monte Carlo simulations of H I line formation in massive YSO 
disk winds 
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ABSTRACT 

Massive young stellar objects (YSOs) are powerful infrared H I line emitters. It has been 
suggested that these lines form in a outflow from a disk surrounding the YSO. Here, new 
two-dimensional Monte Carlo radiative transfer calculations are described which test this hy- 
pothesis. Infrared spectra are synthesised for a YSO disk wind model based on earlier hy- 
drodynamical calculations. The model spectra are in qualitative agreement with the observed 
spectra from massive YSOs, and therefore provide support for a disk wind explanation for 
the H I lines. However, there are some significant differences: the models tend to overpredict 
the "BrafBr-f ratio of equivalent-widths and produce line profiles which are slightly too broad 
and, in contrast to typical observations, are double-peaked. The interpretation of these dif- 
ferences within the context of the disk wind picture and suggestions for their resolution via 
modifications to the assumed disk and outflow structure are discussed. 
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1 INTRODUCTION 

Understanding the formation and pre-main-sequence evolution of 
massive stars (M > 10 Mq) is an important goal in modern as- 
trophysics. Currently, the formation of low-mass stars via gravita- 
tional collapse of a molecular cloud and subsequent disk accretion 
is, at least conceptually, well understood but it is not clear whether 
massive stars form through a scaled-up version of the same process 
or whether more complex, environmental effects are important (see 
e.g. Bonnell, Vine & Bate 2004 for a recent discussion). 

Since massive stars form inside dense molecular clouds, they 
are observed primarily at infrared (IR) and longer wavelengths. 
Embedded massive YSOs are powerful sources of IR H I line emis- 
sion (Simon et al. 1981, 1983; Drew, Bunn & Hoare 1993; Bunn, 
Hoare & Drew 1995; Blum et al. 2004) with equivalent widths 
(EWs) in Bra of ~ 5 - 90 A (Bunn et al. 1995). The lines pro- 
files are complex, displaying both fairly narrow (full-width-at-half- 
maximum [FWHM] ~ 50 - 100 km s _1 ) line cores and broad 
wings extending out to ~ 400 km s _1 in extreme cases. In the 
sample of objects observed by Bunn et al. (1995), the profiles were 
always single-peaked but recent observations reported by Blum et 
al. (2004) include at least one object (NGC 3576) which shows 
clear evidence of double-peaked profiles. By examining flux ra- 
tios across the profiles of lines with differing opacities, Bunn et al. 
(1995) found evidence that, at least in some sources, the lines form 
in an accelerating outflow. 
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Normal, main-sequence OB stars also show H I emission lines. 
These are produced in the star's fast (~ 1000 km s" 1 ) spherical 
wind, but the strength of emission in YSO spectra is much greater 
than in O-stars. If the features in YSOs were due to reprocessing 
of radiation in a spherically symmetric wind, the implied mass-loss 
rates would be up to 10 -6 Mq yr _1 (Simon et al. 1983), substan- 
tially exceeding those of field stars of comparable spectral type. 
Therefore, to interpret the H I observations it is worthwhile to con- 
sider alternatives to formation in a spherically symmetric structure. 

Important observational constraints on the geometry of emit- 
ting material in YSOs are provided by the first overtone bands of 
COat2.3Atm(e.g. Carr 1989; Carretal. 1993; Chandler et al. 1993; 
Chandler, Carlstrom & Scoville 1995). In particular, both Carr et al. 
(1993) and Chandler et al. (1995) conclude that models based on an 
accretion disk surrounding the central object reproduce the CO ob- 
servations of YSOs. Furthermore, although there are some cases in 
which CO data can be modelled in terms of a stellar wind, in gen- 
eral the disk model encounters fewer difficulties (see Chandler et 
al. 1995). The observations strongly suggest that the CO emission 
originates close to the central star: Bik & Thi (2004) and Blum et al. 
(2004) have recently studied CO emission from a range of massive 
YSOs and derived disk radii in the range 0.1-5 AU. Thus it seems 
probable that YSOs harbour accretion disks containing significant 
amounts of hot gas fairly close to the central object and it is natural 
to consider whether the H I emission is associated with such a disk 
or its environment. 

Hamann & Simon (1986) suggested that powerful recombi- 
nation line emission of the sort discussed above originates in an 
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outflow associated not with a stellar wind but with mass-loss from 
a disk. Their work focused on MWC 349, an object whose evolu- 
tionary status is unclear but which possesses spectroscopic signa- 
tures (including H I and CO emission features) broadly similar to 
the embedded massive YSOs (see Hamann & Simon 1986; Kraus 
et al. 2000). Hollenbach et al. (1994) developed a model for out- 
flow from young massive stars by photoevaporation of the outer 
regions of their disks. They were motivated by the need to explain 
the high frequency of occurrence of ultracompact H II regions - 
photoevaporation can provide a source of mass input for H II re- 
gions to balance the loss by pressure-driven expansion. Although 
the rate of mass-loss in photoevaporation models is large, it occurs 
at large radii (typically > 100 AU) and the characteristic flow ve- 
locities are low (comparable to the typical sound speed of 10 - 50 
km s _1 , Hollenbach et al. 1994). Thus, such a flow does not pro- 
vide a promising origin for relatively high velocity features, such 
as the IR H I broad line wings (Bunn et al. 1995). 

Recently, Drew, Proga & Stone (1998) proposed that the in- 
tense radiation field produced by a massive YSO may drive mass- 
loss from the surface of the inner parts of its disk. They performed a 
hydrodynamical simulation which showed that, in addition to driv- 
ing a normal hot star wind component, radiation pressure from a 
massive YSO could propel a dense equatorial flow from a sur- 
rounding disk with terminal velocity in the same regime as the 
observed H I linewidths. These results led them to speculate that 
such a model accounts for the IR H I lines. However, they did not 
perform the radiative transfer calculations required to confirm this 
conjecture. 

To investigate the possibility that massive YSOs harbour disk 
winds which give rise to the observed IR H I line emission, this 
paper presents the results of new two-dimensional radiative trans- 
fer calculations which account for both the complex geometry of a 
YSO disk wind (starting from, but not limited to, the Drew et al. 
1998 model) and the detailed atomic physics of H I line formation. 

In Section 2, the YSO model adopted in this investigation 
and its parameters are discussed. The radiative transfer calculations 
have been performed using the Monte Carlo (MC) code described 
by Long & Knigge (2002) after incorporating a sophisticated treat- 
ment of H I line formation using the approach described by Lucy 
(2002,2003); the method of calculation and code used are discussed 
in Section 3. Results are presented in Section 4 and our conclusions 
discussed in Section 5. 



2 MODEL 

For the radiative transfer calculations (Section 4) a simply parame- 
terised model for the YSO and its disk is adopted. In choosing the 
parameters for this model, our starting point has been to construct 
a reasonable representation of the disk wind obtained in the hydro- 
dynamical simulations presented by Drew et al. (1998). Therefore, 
before discussing our model in detail, a brief description of this 
hydrodynamical model is given. 



of their computational domain. Noting that, for reasonable accre- 
tion rates (~ 10~ 6 M yr _1 ), the luminosity of the central star 
far exceeds the accretion luminosity of the disk, they accounted 
for reprocessing of stellar light by the disk and then computed the 
radiation force using the Castor, Abbott & Klein (1975) parameter- 
isation of the force due to spectral lines. 

Their simulation exhibited a steady state solution with two 
fairly distinct outflow components: a fast (~ 2000 km s _1 ) po- 
lar wind from the central star, rather similar to a normal hot 
star wind; and a slower (< 400 km s _1 ), denser equatorial out- 
flow from the disk at colatitude 6 ~ 60° with mass-loss rate 
~ 3 x 10~ 8 M yr _1 . Between these two outflow components they 
found a complex transitional zone in which the stellar wind stream- 
lines were compressed due to non-radial radiation force compo- 
nents and the presence of the disk wind. 

Drew et al. (1998) suggested that the high density and mod- 
erate velocity of the equatorial disk-wind component make it a 
promising site for the formation of the H I lines observed in massive 
YSO spectra - our goal is to present radiative transfer calculations 
to examine this hypothesis. 



2.2 Geometry 

Figure 1 illustrates the geometrical construction used to describe 
the wind, parameters for which are discussed in the next sub- 
section. In order to capture the essence of the Drew et al. (1998) 
hydrodynamical disk wind model (described above), the structure 
consists of the following four basic elements: 

(i) A central star. The star is assumed to be spherical with radius 

r,. 

(ii) An accretion disk. The disk lies in the :ry-plane of the 
adopted coordinate system and is assumed to extend from the sur- 
face of the star to an outer radius rjisk- For simplicity, it is assumed 
to be geometrically thin and flat. 

(iii) A disk wind. The wind, which is launched from the disk, 
is described following Long & Knigge (2003) and Knigge, Woods 
& Drew (1995): namely the streamlines of the disk wind are as- 
sumed to converge at a point which lies on the 2-axis a distance d 
below the coordinate-origin. Thus the angular boundaries of the 
disk wind are determined by the choice of the inner and outer 
disk radii and the distance d; namely 6 m i„ = tan _1 (r»/d) and 
#max = tan _1 (rdisk/d)(see Figure 1). Motivated by the Drew et al. 
(1998) simulation, the flow is assumed to be stationary and smooth. 

(iv) A spherical stellar wind component. This is assumed to oc- 
cupy the entire polar region above the disk wind and to have radial 
streamlines. The boundary between the stellar wind and disk wind 
is not treated in detail; it is assumed to be perfectly sharp - thus the 
transition zone between the disk and stellar winds which exists in 
the Drew et al. (1998) simulation is neglected here. 



2.1 Summary of the hydrodynamical disk wind model 

Following the methods of Proga et al. (1998), Drew et al. (1998) 
performed a hydrodynamical simulation of a radiatively driven 
wind from a disk around a massive YSO. They considered an early- 
B type star (mass M* = 10 M© , luminosity L* = 8500 L Q , radius 
r, = 5.5 Rq) surrounded by an accretion disk extending from the 
stellar surface to an outer radius rjisk = 10 r», the outer boundary 



2.3 Model parameters 

Table 1 gives a list of the wind parameters adopted for the refer- 
ence model (hereafter, Model A) which are discussed individually 
below. Of our models, Model A and the closely related Model B 
(see Section 4.3) are those most closely matching the Drew et al. 
(1998) simulation. 
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Figure 1. The elements of the geometrical construction used to define the 
wind (only the positive xz-plane is shown - the wind is symmetric about 
both the xj/-plane and the z-axis). The quarter-circle line around the origin 
represents the stellar surface; the stellar wind is shaded light grey, the disk 
wind dark grey and the disk black. The equatorial region beyond the disk 
(white) is assumed to be empty. The figure is not to scale for the parameters 
given in Table 1 . 



2.3.1 Central star 

The parameters of the central star are those of a B1-B2 main se- 
quence star. The effective temperature was computed from the lu- 
minosity using the Stefan-Boltzmann equation. The star is assumed 
to emit radiation as a black body. 



2.3.2 Accretion disk 

In our reference model, raisk = 10 r* is adopted - the same com- 
putational domain considered by Drew et al. (1998). Realistically, 
the disks around massive YSOs are likely to extend to significantly 
larger radii. Therefore, we will also consider disks extending to 
?*disk = 100 r* in Section 4. To consider even larger radii becomes 
impractical since increasing the physical size of the calculation lim- 
its the spatial resolution available in the inner regions - computing 
feasibility prevents us from having arbitrarily many grid cells. It is 
noted that, at least for the treatment of disk radiation adopted in our 
reference model, considering an even larger outer radius will have 
negligible effect - the disk temperature at such large radii will be 
too low for its thermal radiation to contribute significantly to the IR 
spectral range under consideration. 

Given that it lies considerably beyond the scope of this in- 
vestigation to study radiative transfer within the accretion disk in 
detail, approximations must be made regarding the treatment of the 
disk and its emission. In the calculations presented in Section 4, 
results using two different approaches are considered. It is known 
from observations (Henning, Pfau & Altenhoff 1990) that there are 
significant differences in the IR SEDs of massive YSOs - some 
have large near-IR excesses while other have little or no excess. 
Therefore, we have chosen our two models for disk emission to 
approximately bracket the range of plausible IR SEDs for massive 
YSOs. 

In some cases, referred to as a "reflecting" disk, it is assumed 
that the disk effective temperature is determined as a function of ra- 
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Figure 2. The disk temperature (T^sk) as a function of radius (r) for a re- 
flecting disk (solid line) and a reprocessing disk (dashed line). The reflect- 
ing disk temperatures are always lower since they only describe the accre- 
tion luminosity while the reprocessing disk also accounts for re-radiation of 
energy absorbed from the star. 



dius by the accretion rate (exactly as discussed by Long & Knigge 
2002) and that annuli of the disk emit as black bodies at the local 
temperature. During the MC simulations, it is assumed that any ra- 
diation which strikes the disk from above is reflected by a hot disk 
atmosphere rather than absorbed and reprocessed. Since the com- 
putational domain is symmetric about the disk plane, this is closely 
equivalent to the assumption that the disk is completely optically 
thin (transparent). In Model A, the disk is assumed to be "reflect- 
ing". 

Alternatively, as a simple "reprocessing" disk model, it is as- 
sumed that stellar radiation falling on the disk is absorbed and the 
energy re-radiated. In this case, the disk is divided into concentric 
annuli and the energy falling on each annulus from the star com- 
puted following the discussion of Proga et al. (1999). By assuming 
that all this energy is re-radiated locally, and that the annuli radiate 
as black bodies, the disk temperatures are then obtained. The accre- 
tion luminosity is still included in this model but is small relative to 
the reprocessed luminosity. At the start of the MC simulation, the 
disk annuli emit black-body radiation according to their effective 
temperatures - thereby accounting for the luminosity due to repro- 
cessing of stellar light - and during the subsequent propagation of 
MC quanta, any radiation which strikes the disk is simply absorbed. 

Computed disk temperature for both reflecting and reprocess- 
ing disks are shown, as a function of radius, in Figure 2. Spectral 
energy distributions (SED) for radiant energy emitted from mod- 
els with both reflecting and reprocessing disks are shown in Fig- 
ure 3. Comparing the SEDs for reflecting and reprocessing disks 
with rdisk = 10 r» (left-hand panels in Figure 3) shows that the 
treatment of the disk makes a significant difference to the contin- 
uum brightness in the IR region of interest. If the disk only reflects 
(or transmits) light which strikes it, the ~ 2 - 4 /im continuum 
is a combination of stellar radiation and accretion luminosity re- 
leased in the disk. (Figure 3, upper-left panel). However, if the disk 
thermally reprocesses all the stellar light which strikes it, the re- 
processed light is dominant (lower-left panel). Increasing the outer 
radius of the disk to 100 r* increases its brightness. With a reflect- 
ing disk, this effect is fairly small below about 4 /im (upper-right 
panel) because the temperatures of the outer parts of the disk are 
very low. However, since the reprocessing disk has higher tempera- 
tures (Figure 2), changing its outer radius has a greater effect in the 
waveband of interest (Figure 3, lower-right panel). 

We note that since the treatment of the radiation force is grey 
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Figure 3. Spectral energy distributions (SED) for models with star and reflecting disk (upper panels) and star and reprocessing disk (lower panels). The left 
panels are for disks extending to outer radius rj;^ = 10 r* while the right panels are appropriate for r^sk = 100 r*. In each panel, the heavy line shows the 
complete SED while the dashed and dotted lines show, respectively, the contributions due to the star and disk. The same relative flux scale is used in each plot. 
The letters given in the upper right corners of the panels indicate in which of the models discussed in Section 4 the SED is adopted. 



in the Drew et al. (1998) simulation and that in both our approaches 
to the disk SED all the radiative flux which falls on the disk is re- 
emitted locally, both our disk SED models are equally valid in the 
context of the hydrodynamical calculations. 



2.3.3 Disk wind 

The geometry of the disk wind is specified by the distance between 
the origin and the focus point, d (see Figure 1). Here, d = 0.14 r„ 
is adopted leading to # ra j n = 82° and # raax = 89° (see above). The 
dense disk wind component present in the hydrodynamic model of 
Drew et al. (1998) occupied a rather wider angular range (~ 25°). 
However, Proga et al. (1999) showed that improved treatment of the 
radiative line force leads to significantly more swept back equato- 
rial disk winds. In particular, for a model in which the luminosity of 
the central object substantially exceeds that of the disk (their model 
E) they found a disk wind with opening angle 8°. Thus the d-value 
adopted here is chosen to reflect the narrow disk winds obtained by 
Proga etal. (1999). 

The density and velocity structure of the disk wind closely 
follows that adopted by Long & Knigge (2002). Following equa- 
tion (7) of Long & Knigge (2002), it is assumed that the mass loss 
rate per unit area from the disk can be described by 



SM 4q 
oc T dlsk (r) 



(1) 



On the assumption that the mass-loading is proportional to the local 
luminous flux, a = 1 is adopted here. The total mass-loss rate in 
the disk wind is fixed at the mass-loss rate obtained by Drew et al. 
(1998): M = 3 x 10~ 8 M yr -1 . Note that the improved treatment 
of the line force (i.e. Proga et al. 1999) did not significantly change 
the total mass-loss rates from those of the earlier treatment. The 
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Figure 4. The disk wind terminal velocity (Uoo) versus 9 = tan 1 (r/d) 
where r is the disk radius at the base of the streamline. 



poloidal velocity on a streamline in the disk wind is described using 
equation (8) of Long & Knigge (2002), namely 



v(l) 



+ (V 



Rd 

l + Ro 



0d 



(2) 



where c 3 is the sound speed on the disk surface at the base of the 
streamline, Voo is the terminal velocity of the streamline, / is the 
poloidal distance, Rd is a velocity-law scale length for the disk 
wind and /3d is an exponent which determines the acceleration of 
the wind. For simplicity, Voo = w esc is adopted where v e&c is the 
escape velocity on the disk surface at the base of the streamline. 
Figure 4 shows Voo as a function of 6 - these values compare 
favourably with the terminal velocities found by Drew et al. (1998, 
their figure 3). To describe the acceleration of the wind, a velocity 
scale length for the disk-wind, R d = r, and exponent (3 d = 1.5 
are used. 
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Table 1. Parameters for the reference model (Model A). Parameters for 
which departures from the reference values will be discussed are indicated 
with a tick (•/) in the third column. 



Parameter 


Value 


Varied 


Central star 






mass JV/* 


10 Mq 




radius, r* 


5.5 R 




luminosity, L„ 


8500 L 




temperature, T e!! 


2.37 x 10 4 K 




Accretion disk 






inner radius, r m 


r* 


■/ 


outer radius, 


10 r* 


S 


accretion rate, M acc 


10~ 6 M yr" 1 




Disk wind 






focus point, d 


0.14 r* 


S 


mass-loss rate, M 


3 X 10~ 8 M Q yr^ 1 




mass-loading exponent, a 


1 




terminal velocity, Voo 


^esc 




acceleration length, Rd 


r* 




acceleration exponent, f) d 


1.5 




Stellar wind 






mass-loss rate, M s 


10~ 8 M yr- 1 




terminal velocity, v x 


2000 km s" 1 




launch velocity, v c 


10- 2 Doo 




acceleration exponent, f3 s 


1 





of the disk SED has already been discussed (see Section 2.3.2 and 
Figure 3). 

Finally, a model in which an inner hole is introduced to the 
accretion disk will be considered. In addition to increasing r m for 
this model, the d-value is changed in order to preserve the same 
opening angles (#„,;„ and # max ) for the disk wind. 

Departures of the other parameters from their reference values 
will not be considered here. For this investigation we wish to re- 
strict ourselves to central objects with parameters suitable for early- 
type near-main-sequence stars - thus we do not consider models 
with large values of r*. An increase in the stellar radius, r* would 
mimic many of the effects associated with the introduction of an in- 
ner hole to the disk since both push the wind out to regions where 
the rotational velocity is lower. Thus, our model with an inner hole 
can be regarded as a proxy for some of the most important conse- 
quences of an enlarged central star. 

An increase in the stellar luminosity would clearly affect the 
SED, but could also change the mass-loss rate and possibly the 
wind geometry - further hydrodynamical calculations would be re- 
quired to fully investigate results for a range of luminosities. The 
mass-accretion rate (M acc ) is relevant only to calculations with re- 
flecting disks - in such cases, M acc plays an important role in deter- 
mining the disk temperature and IR SED. For reprocessing disks, 
these are controlled by the stellar radiation field and M acc plays 
only a minor role. It is very difficult to reliably obtain M acc by ob- 
servation or theory and so we restrict ourselves to considering only 
one value, that used by Drew et al. (1998). 



2.3.4 Stellar wind 

The density of the stellar wind component is chosen such that it 
would have mass-loss rate M s = 1CP 8 Mq yr -1 if it were spheri- 
cally symmetric. Following e.g. Lucy & Abbott (1993), the velocity 
in this component is assumed to be radial and given by 

v = v c + { Voo - v c )(l - (r./r))" 3 (3) 

Based on the Drew et al. (1998) simulations, = 2000 km s _1 
is adopted. In addition, v c — 10~ 2 t>oo is assumed and /3 S = 1 - a 
standard value in modelling hot star winds - is chosen. 



2.3.5 Departures from the reference model parameters 

In Section 4, the dependence of the results on several of the im- 
portant model parameters listed in Table 1 will be investigated. The 
calculation is, of course, particularly sensitive to the choice of den- 
sity in the wind. The density in turn is sensitive to many of the 
model parameters - most fundamentally, the mass-loss rates which 
determine the mass-loading of the stellar and disk winds. However, 
the velocity-law parameters (e.g. Voo, Rd and (3d) also play a sig- 
nificant role, as does the geometry of the disk wind (determined by 
d). In Section 4, models with density higher than that in Model A 
will be considered: these are created by increasing the mass-loss 
rate only - physically this is the parameter which leads to a simple 
increase in density everywhere in the wind - but it should be noted 
that there is a degree of degeneracy between this and the other pa- 
rameters mentioned above. 

The second departure from the reference model (Model A) 
that will be discussed in Section 4 is an increase in the outer disk 
radius, rjisk- The relevance of this parameter to the determination 



2.4 Atomic data 

For simplicity, it is assumed that both the disk wind and stellar wind 
consist entirely of hydrogen. An atomic model with energy levels 
for principle quantum number n — 1 to 20 and the continuum state 
is used in the calculations. Bound-bound oscillator strengths are 
taken from Menzel & Pekeris (1935). Bound-bound collision rates 
are derived from the oscillator strengths using the van Regemorter 
(1962) formula. Photoionization rates are taken from TOPBASE 
(Cunto et al. 1993) and collisional ionization rates are computed 
using equation 5-79 from Mihalas (1978). In addition to these pro- 
cesses associated with atomic hydrogen, free-free absorption and 
emission by H II ions and Thompson scattering by free-electrons 
are included in the calculations. All other processes are neglected. 



3 METHOD OF CALCULATION 

In the disk wind model, the IR H I line formation is expected to be 
driven by recombination of ionised hydrogen gas fairly close to the 
central star. The observed line ratios (e.g. Bunn el al. 1995) indi- 
cate, however, that the line emission is not optically thin making it 
necessary to perform detailed non-LTE radiative transfer calcula- 
tions to reliably model the line formation process. 

The radiative transfer calculations discussed in the next sec- 
tion were performed using a modified version of the MC code writ- 
ten by Long & Knigge (2002). Since the code has been described 
in detail elsewhere, only a brief overview and discussion of the im- 
portant modifications made for this investigation are presented here 
(see Long & Knigge 2002 for further details of the code). 

To synthesise spectra, the code performs a sequence of MC ra- 
diative transfer calculations. The code uses indivisible energy pack- 
ets as the elementary MC quanta, assumes radiative, statistical and 
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thermal equilibrium in the wind, and utilises a Sobolev treatment 
of bound-bound transitions. 

The code presented by Long & Knigge (2002) adopted a two- 
level approximation in the treatment of line scattering and used ap- 
proximate excitation and ionization formulae for the computation 
of level emissivities. These approximations made the code efficient 
and thus able to handle a large set of atomic data for many chemical 
elements. However, the two-level approximation also meant that it 
was not suitable for modelling lines formed by non-resonance scat- 
tering or recombination. 

Since the IR H I lines form primarily by recombination, the 
code has been substantially modified to incorporate recently devel- 
oped MC radiative transfer techniques which allow the formation 
of such lines to be modelled. In particular, Macro Atoms, as de- 
vised and tested by Lucy (2002,2003), are used. This approach al- 
lows the radiative equilibrium constraint to be rigorously enforced 
at all times without approximation in the treatment of interactions 
between radiation and matter. 

Initially, several MC simulations are performed to compute 
the temperature and degrees of excitation and ionization throughout 
the wind. For these "ionization cycles" the MC quanta are launched 
with frequencies spanning a wide enough range to simulate all the 
important radiative energy inputs to the system. During each of 
these MC simulations, estimators for the radiative heating rate (e.g. 
due to photoionization) are recorded in each grid cell. At the end of 
the simulation, these rates are balanced against cooling rates to de- 
termine the local electron temperature which is then adopted in the 
next simulation. This process is repeated until essentially all grid 
cells pass a chosen convergence threshold in electron temperature. 

For this work, the code has been modified to also record high 
precision MC estimators for the individual radiative rates (pho- 
toionization and bound-bound excitation) following Lucy (2003). 
At the end of each ionization cycle, these are used, together with 
the various collision rates and radiative decay rates, to compute 
level populations and the ionization fraction in each grid cell. These 
populations are used in the next iteration, alleviating the need to 
use the approximate analytic formulae for ionization and excita- 
tion employed by Long & Knigge (2002). No convergence criterion 
is enforced for the level populations at present. Lucy (2002) has 
shown that, provided the populations of the lowest levels are esti- 
mated sensibly, the indivisible-energy-packet MC method produces 
highly accurate emissivities even if the populations for the emitting 
upper levels are relatively poorly known. This insensitivity moti- 
vates our use of the MC method since it substantially reduces our 
reliance on detailed calculation of excited level populations for an 
adequate non-LTE investigation of recombination line formation. 
We note that our convergence criterion on the electron temperature 
implies a corresponding convergence in the energy flow between 
the radiation field and the thermal energy pool. Since this energy 
flow is a summation of all the radiative heating processes in the hy- 
drogen model, its behaviour implies good convergence amongst all 
the processes which are energetically important. 

After the "ionization cycles", a set of additional MC simu- 
lations are performed to compute the spectrum (the "spectral cy- 
cles"). During these cycles the MC quanta are created in only a 
fixed frequency interval while the thermal, ionization and excita- 
tion states of each grid cell are fixed to the values previously deter- 
mined as described above. By tracking the quanta, we compute the 
spectrum as seen by observers at different inclination angles in the 
frequency range of interest. 

In all the calculations described below, the model structure 
and radiation field properties are discretized onto a two-dimension 



Table 2. The observed range of properties of massive YSO H I IR lines 
in the sample discussed by Bunn et al. (1995). Note that in several of the 
objects, Br7 and/or Pf7 are not observed and so the range of EW for these 
lines is biased towards the brighter sources compared with Bra. 



Line 


-EW 


FWHM 


HWZI 




(A) 


km s -1 


km s~ x 


Bra 


4.8 - 87 


55-155 


50-415 


Br7 


2.4 - 23 


70-150 


70 - 270 


Pfy 


2.3-12 


110-150 


40 - 260 




EWR 




EWR 


Bra/Br7 


1.1-4.4 


Bra/Pf7 


3.0-7.2 



(r-z) grid consisting of 3000 cells. During each MC calculation, 
the input radiative energy (which arises from the star and accretion 
disk) is divided into 3 x 10 6 indivisible quanta. 

We have explicitly verified that our approach is adequate for 
the calculation of the spectrum by comparing results for models 
with identical input parameters using smaller and larger numbers of 
grid cells. These calculations differed by at most a few per cent in 
the computed H I line strengths from the standard (3000 cell) case. 
This implies that the number of packets, number of grid cells and 
the quality of the MC estimators are all sufficient for calculations 
to this degree of accuracy. 



4 RESULTS 

4.1 Observational constraints 

Early observations of IR H I emission lines in massive YSOs (e.g. 
Simon et al. 1981) had insufficient signal-to-noise ratio to provide 
detailed information on the line shapes, thus subsequent theoretical 
work (e.g. Hoflich & Wehrse 1987) was concerned primarily with 
modelling the observed line fluxes and flux ratios. 

More recently, higher quality data (e.g. Drew et al. 1993, Bunn 
et al. 1995) has allowed the line shapes to be studied in greater de- 
tail in several massive YSOs. In common with earlier work (e.g. 
Simon et al. 1981,1983; Drew et al. 1993), Bunn et al. (1995) ob- 
served a moderately strong, pure-emission line of Bra in all of their 
targets. In most cases, they also reported weaker emission in Br7 
and Pf7. Their observed line profiles are single peaked and not 
shifted from the expected local rest velocity by more than a few 
km s _1 . Many of the profiles show line wings which extend to sev- 
eral hundred km s~ 1 from line centre and in some cases these wings 
are clearly asymmetric. By considering wavelength-dependent line 
ratios (following Drew et al. 1993), Bunn et al. (1995) found evi- 
dence, in at least some cases, of hybrid profiles consisting of a nar- 
row line core (which may have a nebular origin) and broad wings in 
which the line ratio is consistent with formation in an accelerating 
outflow. 

High resolution observations of several massive YSOs have 
very recently been presented by Blum et al. (2004). Interestingly, 
in their sample, at least one object (NGC 3576) has double-peaked, 
narrow line profiles. However, it is difficult to study the broad line 
wings in the Blum et al. (2004) data since the spectral coverage is 
too narrow to establish the underlying continuum with certainty. 

In the following sub-sections, results of our models will be 
compared primarily with the observations presented by Bunn et 
al. (1995) since their spectra provide good constraints on both line 
strengths and shapes for a significant sample of objects. Unfortu- 
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Figure 5. Model A (reference model with reflecting disk): computed Bra (upper), Br7 (centre) and Pf7 (lower) line profiles for viewing angles to the polar 
axis (from left to right) of 30°, 45°, 60°, 70°, 80° and 85°. The flux is given for an unreddened source at a distance of 1 kpc and the velocity is measured 
relative to line centre. 



nately, owing to the great difficulty in establishing the extinction 
along the line of sight to massive YSOs, interpreting line fluxes 
and flux ratios for lines at different wavelengths is very difficult (it 
is known from Bra/Pf/y flux ratios that the H I lines are not well 
described by either the LTE or the Baker Menzel case B recombi- 
nation assumptions - thus the Bra/Bry ratio cannot be trusted as an 
indicator of the reddening [Hoflich & Wehrse 1987]). Therefore, in 
our discussion, we prefer to focus on EWs, EW ratios (EWRs) and 
line shapes since these are not affected by reddening. For quantita- 
tive comparison, Table 2 gives a summary of the range of observed 
EW, FWHM and half-width-at-zero-intensity (HWZI) from Bunn 
et al. (1995). The table also gives the observed range of EWRs for 
Bra/Bry and Bra/Pf/y 1 . Since our model parameters are not fine- 
tuned to one particular object, these values should be regarded only 
as indicative of the appropriate regime in which acceptable model 
predictions should lie. 

4.2 Reference model (Model A) 

Spectra have been computed in the 2-5 fim range using the refer- 
ence wind model (Model A) which was described in Section 2. In 
this wavelength range, the model predicts moderately strong emis- 
sion in the first three lines of the Brackett series and weaker emis- 
sion in lines of the Pfund series (Pfy3 and Pfy). This is in qualitative 
agreement with the observed IR properties of massive YSOs (e.g. 
Simon et al. 1981,1983; Drew et al. 1993; Bunn et al. 1995). 

1 Hereafter, these EWRs will be referred to as EWR[Bra/Br7] and 
EWR[Bra/Pf7] respectively. 



Figure 5 shows the predicted Model A line profiles of Bra, 
Br7 and Pf/y. Spectra are shown for a grid of viewing angles rang- 
ing from 30° to 85° (measured relative to the polar axis). For 
quantitative comparison with observations, the computed EW and 
FWHM of these lines are tabulated in Table 3. For completeness, 
computed line fluxes and flux ratios are given in Table 4. 

For all three lines, the EW is largest for viewing angles (0 O bs) 
of 70 - 80° and becomes noticeably smaller at lower O bs- The com- 
puted Bra EWs are consistent with the high end of the range of 
measured EW (Table 2) obtained by Bunn et al. (1995). Similar 
agreement between observations and the model is also found for 
Pf/y making the modelled EWR[Bra/Pfy] similar to the observed 
ratio. 

For Bry, the Model A computed EWs lie in the mid-range 
of the observed values (see Table 2). The relative weakness of the 
computed Br7 line means that the model EWR[Bra/Br7] is sys- 
tematically too high - in the observations the EWR for these lines 
ranges from 1. 1 to 4.4 but the computed Model A ratio is somewhat 
higher, ranging from 11 to 13 for the viewing angles considered. It 
is probable that part of this disagreement is the result of opacities 
being too low in the models - but given the good agreement with 
observations for EWR[Bra/Pf/y], it seems unlikely that a simple 
underestimate of the opacity in Bra is solely responsible. Perhaps 
more importantly, the EWR[Bra/Bry] will be affected by inaccu- 
racies in the model continuum shape - since Bra and Pf/y lie at 
similar wavelengths, the choice of continuum shape will have little 
effect on the computed ratio of their EWs but because Bry lies at a 
significantly shorter wavelength, the EWR[Bra/Br7] is much more 
sensitive. An alternative treatment of the continuum - that invoking 
a reprocessing disk - is considered in Section 4.2. It is also noted 
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Table 3. Equivalent widths (EW) and full-width-at-half-maximum (FWHM) for the Bra, Br7 and Pf/y lines at different viewing angles (9 ohi ). The values of 
FWHM have accuracy of around ~ 50 km s _1 , limited by the frequency gridding of the calculations. The accuracy of the EWs is limited by MC noise in the 
spectrum. Typically, the Monte Carlo noise in the continuum is < 2 per cent. 
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Model C (as Model B but with mass loss rate enhancement) 
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Model D (as Model A but with disk out to 100 r„) 
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Model E (as Model C but with disk out to 100 r*) 
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Model F (disk with inner hole, see Section 4.4) 
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a For Model B, the Br 7 and Pf 7 lines are very weak and only appear marginally above the MC noise - therefore, the entries for these lines give only upper 

limits on their EW. 



that, compared with Bra and Pf 7 , observations of Br 7 for heavily 
embedded targets are more likely to be contaminated by emission 
from other sources (e.g. larger scale nebular emission) owing to the 
wavelength dependence of the source extinction. 

The computed line profiles are complex, usually exhibiting 
double-peaks which become more pronounced at larger (9 bs- This 
is in contrast to the observations where only single-peaked pro- 
files are usually observed (see Section 4.1). Also, the computed 
profiles are too broad (FWHM ranging from 300 to 600 km s _1 
as a function of # b s ). P art °f this discrepancy between model and 
observations may be explained by contamination of the observed 
profiles by narrow nebular emission - there is some evidence for 
hybrid line profiles with narrow cores (see Section 4.1) - however, 
the half-width-zero-intensity (HWZI) measurements also suggest 
that the observed line profiles do not extend more than 200 to 300 
km s _1 from line-centre in most cases. 

The large line widths and double-peaked profiles arise due to 
the rotational velocities in the region of line formation. Figure 6 
identifies the region of line formation for the Bra, Bry and Pf 7 
lines in Model A; specifically, the plots show the amount of energy 
that escapes in each spectral line from each grid cell in the model. 
Note that MC noise is responsible for the graininess of the figures. 
It can be seen that in all three cases, the region of line formation 
is close to the disk surface, within the first few velocity law scale- 
lengths along the wind streamlines. Bra tends to form slightly fur- 



ther out than either Br 7 or Pf 7 - this is to be expected owing to the 
higher opacity in the Bra line. 

To obtain models which produce narrower lines, either more 
low-velocity gas must be added, or high-velocity gas must be re- 
moved. The simplest way in which to introduce low- velocity gas is 
by increasing the outer disk radius (rdisk) - models with wider disks 
are presented in Section 4.3. To remove high-velocity gas requires 
a reduction in the wind density at small values of r. Given that the 
inner disk is brighter than the outer disk, it seems unlikely that a 
radiatively driven flow will not have greatest mass-loading at small 
radii. Therefore, the most plausible way in which to eliminate high- 
velocity gas is not to change the prescription of mass-loading but 
rather to consider models in which either the stellar radius is larger 
than that adopted in the reference model or the disk inner radius 
is several times greater than the adopted stellar radius. An exam- 
ple of the latter, a model with an empty inner cavity in the disk, is 
discussed in Section 4.5. 

To assess the relative importance of the disk wind and stellar 
wind component, results from a calculation considering only the 
disk wind component were compared with those which include the 
stellar wind. These two calculations give virtually identical results 
- this is not unexpected since it is already known that a spherical 
wind would require a mass-loss rate at least an order of magni- 
tude greater than that in the stellar wind component here (Simon 
et al. 1981, see Section 1). In view of this, in all the calculations 
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Table 4. Computed Bra line fluxes and flux ratios Bra/Br7 and Bra/Pf'7 at different viewing angles. Fluxes are given in units 10 13 ergs cm 2 s 1 for an 
unreddened source at a distance of 1 kpc. 
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1.51 


1.56 


1.64 


Bra/Pf'7 3.07 3.44 3.66 


4.11 


4.57 


4.90 
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a For Model B, the Bry and Pfy lines are very weak and only appear marginally above the MC noise - therefore, the entries for these lines give only lower 

limits on the appropriate flux ratios. 



discussed below, the stellar wind component is dropped for com- 
putational efficiency. 



4.3 Reprocessing of radiation by the disk (Models B and C) 

In Model A (discussed above) it was assumed that the photons 
striking the accretion disk were reflected rather than absorbed. To 
investigate an alternative simple hypothesis, namely that the disk 
absorbs and thermally re-emits the radiation which strikes it (see 
Section 2.2.2), a second model (Model B) is now considered. This 
is the extreme, optically thick disk case. All the system parameters 
(Table 1) for Model B are identical to those of Model A, saving that 
the stellar wind component is omitted (see above). 

The simple treatment of reprocessing adopted (see Section 
2.2.2) means that the disk is much brighter at IR wavelengths in 
Model B than Model A. Thus, although the wind parameters in 
Model B are unchanged, the EWs of all the emission lines are much 
smaller (the EWs are given in Table 3). At moderate to high # obs , 
the Model B EWs for Bra and Pf7 are comparable to the observed 
EWs in the weaker sources reported by Bunn et al. (1995). How- 
ever, the computed Br7 line is very weak. 

To increase the line EWs in the presence of a reprocessing 
disk the wind density (and hence emission measure) needs to be in- 
creased. Therefore, a further model (Model C) has been constructed 
- it is identical to Model B saving that the disk-wind mass-loss rate 
has been increased by a factor of four. A higher mass-loss rate may 
be feasible even under the assumption of driving purely by radiation 



since the absolute mass-loss rates provided by the hydrodynamical 
model of Drew et al. (1998) may be uncertain by a factor of a few 
owing to the parameterisation of the line-driving force. EWs and 
FWHM for Model C are given in Table 3 and the line profiles com- 
puted from this model are shown in Figure 7. 

The Model C EWs are slightly greater than even the largest 
observed values from the Bunn et al. (1995) data, suggesting that 
the factor of four mass-loss rate increase between models B and 
C comfortably brackets the range of values that are likely to be 
consistent with the data for particular objects. 

The EWR[Bra/Pf7] is smaller for Model C than Model A, ly- 
ing close to the mid-range of observed values. The EWR[BrafBry] 
is also smaller in Model C - it is still greater than the observed ratio 
in any of the objects observed by Bunn et al. (1995) but the discrep- 
ancy is much less. This change in the line ratios is mostly due to 
the increased Bra opacity resulting from the higher wind density. 
This result suggests that by further increasing the mass-loss rate 
an EWR[Bra/Br7] close to that observed could be achieved. How- 
ever, it is noted that a further increase in the wind density will lead 
to Bra EWs which are uncomfortably large. 

The Model C line profiles are generally less double-peaked 
and more square-topped than the Model A profiles. Their FWHM 
are very similar to those in Model A and the line wings still extend 
to several hundred km s _1 at the base. 

In general, Model C appears to suggest that the combination 
of a reprocessing disk and a higher disk-wind mass-loss rate is in 
closer agreement with the observations than the reference model 
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A - C by adding low-velocity gas. In this section two more models 
(Models D and E) are considered: these models are identical, re- 
spectively, to Models A and C saving that the outer disk radius is 
set to 100 r* rather than 10 r*. 

The computed line profiles for Model D are shows in Figure 8 
and the EWs for Models D and E are given in Table 3. The dif- 
ferences between the profile shapes shown in Figures 5 and 8 are 
generally fairly small. The double-peaked shape is slightly filled 
in when the larger disk is considered; however, the assumption that 
the mass-loading of streamlines is proportional to the luminous flux 
(see Section 2.2.3) means that the majority of the mass-loss still 
originates from the inner disk where velocities are high. The most 
important consequence of the larger disk is the increase in the IR 
continuum level (see Figure 3 and the discussion in Section 2.3.2). 
The brighter continuum means that, when the disk is bigger, the 
line EWs are all smaller. Also, the EWR[Bra/Br7] is smaller - the 
larger disk radius increases the continuum around 4 /im by more 
than at 2 /im (see Figure 3). 

For Model D, the Bra EWs lie comfortably in the mid-range 
of the observed values and the EWR[Bra/Pf7] is similar to that 
found in Model A. The EWR[Bra/Br7] is smaller in Model D than 
A (as required) but the improvement is less significant than ob- 
tained with Model C (see Section 4.2). Model E predicts smaller 
EWs, but still within the range observed. The EWR[Bra/Br7] in 
Model E is closer to that observed than in any of the other mod- 
els - this follows since Model E benefits from the improvement in 
this ratio resulting from both a higher line opacity (density) and 
stronger disk continuum. 
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Figure 6. Regions of line formation for Bra (upper), Br7 (centre) and Pfy 
(lower). Note the logarithmic axes. The solid black lines define the bound- 
aries of the disk wind (the star surface, the disk surface and the inner and 
outer wind boundaries defined by 8 m [ n and 6> max ). Within the wind, the vol- 
ume is divided into the computational grid cells which are shaded propor- 
tional to the energy escaping to infinity from that cell in the spectral line 
under consideration. The shading scale is logarithmic. The graininess is due 
to MC noise. 



(Model A). But further hydrodynamical modelling is required to 
examine whether so high a mass-loss rate is feasible. 



4.4 The influence of the outer radius of the disk 

The outer disk radius paralleling that of Drew et al. (1998) and 
adopted in Models A - C is, arguably, rather small. In principle, 
line formation in a wind rising from the outer parts of a larger disk 
could eliminate the double-peaked line shape obtained with Models 



4.5 The influence of a central hole in the disk 

It is apparent that a failing of Models A - E in comparison with 
the data is that the computed line profiles are too broad. One of 
two possible solutions to this problem is to invoke an inner hole in 
the accretion disk, thereby removing the highest velocity material 
from the region of line formation. Optically thin inner cavities in 
accretion disks around young luminous stars have previously been 
suggested in various contexts (e.g. Herbig Ae/Be stars - Dullemond 
et al. 2001; Tuthill, Monnier & Danshi 2001; Natta et al. 2001; and 
see Monnier et al. 2005 for a recent discussion). Interestingly, a disk 
with an inner hole might also help explain the residual discrepancy 
between the observed EWR[Bra/Br7] and the models discussed 
above. Since the inner parts of the accretion disk are hottest, re- 
moving them will push the SED of the disk to the red. This will 
increase the continuum level at Bra relative to Br7 and therefore 
reduce the EWR[Bro/Br7] as required. However, a quantitative in- 
vestigation of this possibility goes beyond the scope of this paper 
since more sophisticated models for the disk emission than used 
here would be required. 

Consideration of models in which the disk has a central hole 
represents a significant departure from the Drew et al. (1998) model 
and therefore, until further hydrodynamical calculations are per- 
formed, it is difficult to place useful constraints on the likely ge- 
ometry and parameters for mass-loss from such a disk. Neverthe- 
less, for illustrative purposes, a model adopting the same wind 
opening angles (# m i n , #max) and mass-loss rate as the reference 
model (Model A; Section 4.1) but with an inner disk radius of 
ri n = 5r, and outer disk radius of r- m — 100 r* has been con- 
structed (Model F) - this corresponds to a disk with a circular inner 
hole of ~ 0.1 AU. It is noted that, conceptually, invoking a wind 
launched significantly further out than the stellar surface is reminis- 
cent of the photoevaporation model (e.g. Hollenbach et al. 1994) - 
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Figure 7. As Figure 5 but showing Model C (reprocessing disk, enhanced mass-loss rate) results. 
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Figure 8. As Figure 5 but showing Model D (reflecting disk with larger outer radius) results. 
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Figure 9. As Figure 5 but showing Model F (disk with inner hole) results. 



however, the distance from the central star at which the mass-loss 
occurs here is still much smaller than expected from photoevapo- 
ration. It is, however, comparable to the smallest disk radii inferred 
from CO observations by Bik & Thi (2004). 

The computed FWHM and EWs for Model F are tabulated in 
Table 3 and the profiles are shown in Figure 9. The EWs are smaller 
than those for Model A because the density in the wind is lower, a 
result of the radial dependence of the volume element in a cylin- 
drical polar system. The lower density also leads to lower opacity 
and hence larger ratios of Bra to Br7 and Pf7. However, as antic- 
ipated, the Model F lines are narrower than in Model A, by up to 
100 km s _1 in FWHM for moderate inclination angles. The profile 
widths at base are also smaller - in Model F the profiles cut-off at 
~ 300 km s _1 from line centre. Thus holes of the size considered 
here, or a little larger, could explain why the HWZI measurements 
of Bunn et al. (1995) are rather smaller than would be expected 
based on Models A - E. 

Given that there are no strong constraints on the hole size, ge- 
ometry or mass-loss rate appropriate for the model here it is not 
possible to drawn definitive conclusions. However, it seems that 
models with an inner disk cavity have the potential to resolve the 
discrepancies between the computed disk wind line profiles and 
the observations. Model F also suggests that the second possibility 
mentioned in Section 4.2 - that of an enlarged stellar radius - can 
also provide a possible solution since increasing the stellar radius 
will affect the line widths in a similar way to introducing an inner 
gap in the disk. Constraints on the size of any inner hole, the stel- 
lar radius and the geometry of the disk wind are needed in order to 
proceed. These could be provided via sophisticated observational 
techniques such as interferometry (which has already been applied 



to Herbig Ae/Be stars by e.g. Monnier et al. 2005) or spectropo- 
larimetry (as discussed by Vink, Harries & Drew 2005). 



5 DISCUSSION 

As mentioned in Section 4, direct confrontation of theoretical line 
fluxes with observations is difficult owing to the uncertainty in the 
extinction along lines-of-sight to massive YSOs. But on comparing 
EWs - the quantities most readily extracted from observations - the 
radiative transfer calculations presented above lead us to conclude 
that model disk winds, such as that obtained by Drew et al. (1998), 
do predict H I line strengths consistent with observations of massive 
YSOs. This matching is achieved for total mass-loss rates between 
0.3 and 1.2 x 10~ 7 M© yr _1 , in contrast with previously calculated 
spherical models for early B stars which require mass loss rates up 
to 10~ 6 M e yr" 1 (Simon et al 1983, Nisini et al 1995). 

There are, however, some difficulties when the models are 
considered in detail. Although the computed EWR[Bra/Pf7] 
agrees well with observations, the EWR[Bra/Br7] tends to be over- 
predicted by the models. We have investigated the effect of a higher 
density on this ratio by considering models with greater mass- 
loss rate (e.g. Model C; recall that the model density is sensitive 
to several parameters in addition to the mass-loss rate and thus 
there is a degree of degeneracy among these quantities - see Sec- 
tion 2.3.5). As expected, increasing the wind density decreases the 
EWR[Bra/Br7] owing to the greater optical depth. However, to 
fully explain the EWR discrepancy solely by such a modification 
is difficult since the absolute values of the EWs grow rapidly if the 
density is raised. At the same time, this EWR is very sensitive to the 
assumed underlying SED, owing to the significant wavelength dif- 
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ference between the Bra and Br7 lines (unlike Bra and Pf7 which 
lie at similar wavelengths). Observationally it is clear that there is 
genuine diversity in the IR SEDs of massive YSOs (e.g. Henning 
et al. 1990). In our models we have considered only two limiting 
cases for the treatment of the SED and have not concerned our- 
selves with trying to fit specific objects. To determine whether the 
EWR[Bra/Br7 discrepancy is unavoidable, there is now a need to 
move on to fitting observations of individual well-observed sources 
with reliable reddening estimates. 

Important differences are also evident in comparisons of the 
line profiles obtained from our models and those typically observed 
in massive YSOs. To illustrate this, Figure 10 compares the ob- 
served Bra profile for GL989 (taken from Bunn et al. 1995) with 
computed profiles from two of the models. GL989 is chosen here 
since the line profiles for this object have strong extended wings 
and, in contrast to several of the other objects in the Bunn et al. 
(1995) sample (e.g. S106IR), have simple profiles suggesting a sin- 
gle dominant region of line formation. The comparison models and 
viewing angles (Model B with 6> b s = 60° and Model E with 
#obs = 45°) were selected since they predict Bra EWs close to 
the observed value, thereby allowing a direct comparison of profile 
shape. 

The figure shows that the Model B Bra profile at # bs = 60° is 
too broad and, in contrast to the observations, is double-peaked. 
The Model E profile (6U S = 45°) is in better agreement with the 
observation - the combination of higher optical depth in the model 
(owing to the increased mass-loss rate) and lower viewing angle 
reduce the profile width and suppress the double-peaks - but a 
significant discrepancy remains. To resolve this difference in line 
shape between the models and typical observations, one could in- 
voke even smaller inclination angles. Although viewing angles are 
not known for most luminous YSOs, this is not an attractive solu- 
tion since small viewing angles (close to pole-on) are statistically 
disfavoured in an unbiased sample. Furthermore, at least one of 
the objects discussed by Bunn et al. (1995) is known to be viewed 
almost edge-on (S106IR; Solf & Carsenty 1982). Thus we prefer 
to pursue resolutions to the two styles of line shape discrepancies 
which do not rely on orientation effects. 

Increasing the physical extent of the disk from which the wind 
is launched does not readily eliminate the double-peaked line shape 
from our models (Models D and E) - although some of the wind 
does occupy regions with lower rotational velocities when the disk 
is made larger, our assumption that the mass-loading is propor- 
tional to the luminous flux of the disk means that most of the mass 
in the wind is still launched at small radii. The problem of recon- 
ciling observed single-peaked line profiles with theoretical predic- 
tions of double peaks has arisen elsewhere, particularly in studies 
of nova-like variables (see e.g. Home 1997) and active galactic nu- 
clei (AGN). In the context of AGN, Murray & Chiang (1997) have 
argued that substantial radial velocity shear in a disk wind can sup- 
press the formation of double-peaked profiles as required. For this 
to occur, it is necessary that line optical depths are high such that 
photon escape is strongly favoured in the poloidal direction where 
accelerating outflow ensures steep velocity gradients. It is this dif- 
ference between Model C, with enhanced mass loss, and Model A 
that explains the relatively flat-topped profiles in Fig. 7, compared 
to Fig. 5. Hence there is scope to re-organise model parameters in 
order to exploit this mechanism to suppress double peaks. In par- 
ticular, a way of raising line opacity without also raising Bra EWs 
is needed. This will warrant further investigation in future work, 
focusing on fitting observations of specific objects. 

Some of the discrepancies, in particular the large linewidths 



predicted by the models, may be explained by departures from the 
particular disk wind geometry adopted (which was motivated by the 
Drew et al. 1998 model). Two interesting possibilities are those of 
an enlarged radius for the central star, relative to a main-sequence 
radius, or of an inner gap in the YSO disk. These possibilities have 
rather similar consequences - by moving the disk wind out to larger 
radial distances from the central object, the correspondingly lower 
rotation velocities reduce the extent of the line wings and make 
the double-peaks less prominent (this is illustrated by our model 
with a cavity in the inner disk, Model F). Also, the SED for a disk 
with a hole will be redder than that for a complete disk, helping to 
explain the EWR[Bra/Br7] discrepancy. These alternatives might 
be distinguished by observational techniques that probe the inner 
disk geometry (e.g. interferometry and spectropolarimetry) while 
the influence of a gap on the SED needs to be investigated fully 
in subsequent work incorporating more sophisticated treatments of 
the disk SED than used here. 

A third possibility is that of disk winds in which the stream- 
lines diverge less at large distances. In our models, the divergence 
of the streamlines is imposed by our adopted "split dipole" geom- 
etry. In a geometry with less divergence of the flow, the density 
would be higher in the outer parts of the wind for a given mass- 
loss rate and velocity law. The increased emission measure at large 
radii would help to wash out the double-peaked profiles which form 
when the emission is dominated by gas in rotation close to the 
central object. For appropriate viewing angles, line-of-sight opti- 
cal depths would also be larger for a less divergent outflow. As has 
already been noted, higher optical depths present a promising so- 
lution to part of the EWR[Bra/Br7] discrepancy. Thus, in the con- 
text of studying particular objects in detail, considering a range of 
different disk wind geometries is likely to prove fruitful in future 
work. 

These issues aside, the major success of our radiative transfer 
simulations is their support of the hypothesis that disk winds can 
play a significant role in creating the H I lines of massive YSOs. In 
particular, the departure from one-dimensional wind models has al- 
leviated the need to invoke the uncomfortably large mass-loss rates 
required in earlier spherical models. Further investigations of the 
applicability of the disk wind model to particular objects are now 
warranted. 
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